Method and apparatus for imaging artefact reduction

ABSTRACT

The invention relates to a method and apparatus for reducing imaging artifacts in physiological images, such as magnetic resonance imaging images. A model of a physical cause for the image artifact is modelled and the effects of perturbing this model on a focus criterion of the image analyzed so as to optimize the focus criterion. For an optimized focus criterion, the model should reflect the actual physical cause of the image artifact. A focus criterion is a measure of entropy of the image, with an optimized image having reduced entropy.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to image processing and in particular toimage processing techniques whereby image artefacts produced duringimage acquisition are corrected for in the generation of an outputimage.

2. Discussion of Prior Art

Magnetic Resonance Imaging (MRI) is a widely used technique for medicaldiagnostic imaging. In a conventional MRI scanner, a patient is placedin an intense static magnetic field which results in the alignment ofthe magnetic moments of nuclei with non-zero spin quantum numbers eitherparallel or anti-parallel to the field direction. Boltzmann distributionof moments between the two orientations results in a net magnetisationalong the field direction. This magnetisation may be manipulated byapplying a radiofrequency (RF) magnetic field at a frequency determinedby the nuclear species under study and the strength of the appliedfield. In almost all cases, the species studied is the nucleus of thehydrogen atom, present in the body mainly in water molecules, and RFpulses are applied at the resonant frequency of these water protons.

The energy absorbed by nuclei from the RF field is subsequentlyre-emitted and may be detected as an oscillating electrical voltage, orfree induction decay (FID) signal in an appropriately tuned antenna.More commonly, a further RF pulse, or magnetic field gradient, is usedto postpone signal acquisition and generate a spin-echo or gradient-echosignal.

Spatial information is encoded into the echo signal by virtue ofadditional linearly varying magnetic fields, known as magnetic fieldgradients, applied during or prior to echo acquisition. The principle ofspatial encoding is that in the presence of the field gradient, the netfield experienced by a given nuclear moment, and hence its resonantfrequency, is a function of position within the scanner. When a gradientis applied during the echo acquisition, the received signal contains arange of frequency components representing nuclei at different locationsalong the gradient direction. Fourier transformation of this signalyields a one-dimensional projection through the patient. This techniqueis known as frequency encoding. Two-dimensional encoding requires use ofan additional gradient applied perpendicular to the frequency encodingaxis, known as the phase encoding gradient. This gradient is applied fora short time prior to data acquisition. The acquisition process isrepeated perhaps 256 or 512 times using phase encoding gradients ofdifferent strengths. Simultaneous frequency and phase encoding yields atwo-dimensional data set which when subjected to two-dimensional Fouriertransformation provides the required image. This array of data exists inwhat is known as k-space and is the Fourier transform of the imagespace. The effect of the phase encoding gradient is to move the start ofthe data acquisition to a particular location along one axis in k-space(dependent on the gradient strength), whilst frequency encodingrepresents a sweep through k-space parallel to the other axis. Each ofthese sweeps is known as a “shot” or “view”.

Spatial localisation in the third dimension may be achieved using anadditional phase encoding gradient, or more commonly by using a gradientand narrowband RF pulse to restrict the initial perturbation of nuclearmoments to a single tomographic slice. This principle can readily beextended to multislice MRI.

In conventional MRI, a single phase-encoding view is acquired after eachRF excitation. However, faster imaging sequences now exist in whichfurther RF pulses and phase encoding gradients are used to acquire atrain of differently encoded echoes after each excitation. These echoestraverse several lines of k-space and reduce scanning time by a factorequal to the echo train length. In the extreme case, single shot echoplanar imaging (EPI) techniques cover the whole of two dimensionalk-space in a single acquisition lasting less than 100 ms, althoughspatial resolution and image quality may be significantly compromised.

Patient movement during the acquisition of MRI images results indegradation of the images that can obscure the clinically relevantinformation. Each readout period takes a few milliseconds (ms), whereasthe time interval between readouts might be between 500 and 2000 ms. Themajority of blurring and ghosting artefacts caused by patient motion aredue to motion between lines in k-space, rather than motion during asingle readout.

Movement leads to phase errors between lines of k-space, which in theresulting image appear as blurring and ghosting along the phase encodedirection. These phase errors can result from translations androtations. Translations of the patient in the readout direction resultin a frequency dependent phase shift in each line of k-space. Rotationsin the spatial domain are also rotations in k-space, and result in aphase shift that is a more complicated function of position.

A particular type of MRI image investigation, known as diffusionweighted imaging, takes place in the presence of an additional andseparate gradient. The integral over time of the diffusion weightedgradient is greater than the integral for the phase encode or readoutgradients. The purpose of this additional gradient is to make the imagessensitive to molecular motion of the order of 10 μm. A side effect isthat the images are also sensitive to bulk motion on the same scale.Anderson et al. in Magn. Reson. Med. Volume 32, 1994, pages 379-387 haveshown that for small rigid body movements, the resulting artefacts canbe modelled using zero and first order phase correction terms.

Considerable work has been done by MR researchers to model patientmotion, and to attempt to correct for it. The impact of different typesof motion on the resulting images is well understood, but clinicallyusable retrospective motion correction techniques are not yet available.Existing algorithms tend to correct only for one dimensional motion, orthey require exotic image acquisition strategies that are not generallyapplicable.

There are broadly two classes of movement correction algorithm used inMRI; with and without “navigator echoes”. Motion correction usingadditional echoes referred to as “navigator echoes” involves theacquisition of additional echoes that are not phase encoded between eachphase encoded echo. All navigator echoes are projections through theobject. It is therefore possible to measure motion between the navigatorechoes, and consequently infer the motion between the correspondingphase encoded echoes. The navigator echoes are most commonly used tomeasure motion in one or more translational directions. Published papersdescribing the use of navigator echo techniques include that of Ehman etal. in Radiology, Volume 173, 1989, pages 255 to 263. Recently severalauthors have proposed obtaining rotational information either from thesestraightforward navigator echoes, for example Anderson et al. in Magn.Reson. Med. Volume 32, 1994, pages 379-387, or from circular navigatorechoes for example Fu et al. in Proc. Soc. Magn. Reson., 1994, page 355.

There have been attempts to measure movement directly from the phaseencoded data. Felmlee et al., as reported in Radiology, Volume 179,1991, pages 139 to 142 tried to measure translations directly from ahybrid space comprising the Fourier transform of the readout vs. phaseencode, but found that it worked for phantoms with high spatialfrequency edges, but on human subjects only if high contrast markerswere used. A possible solution to this is to acquire spiral readouts,all of which sample a range of spatial frequencies as described byKhadem et al. in Proc. Soc. Magn. Reson., 1994, page 346, but this isimpractical on the majority of MRI hardware. An alternative strategydescribed by Wood et al. in J. Magn. Reson. Imag., Volume 5, 1995, pages57 to 64, is to locate discontinuities in k-space that correspond tosudden movements of the patient, to split the regions of k-space betweenthese discontinuities into sub-images, then to correct for translationsby applying a phase shift to k-space, and rotations by rotating in thespatial domain. This technique cannot correct for continuous movementduring acquisition, and would appear to require considerable userinteraction.

An alternative approach to motion correction without navigator echoes isthe Projection Onto Complex Sets (POCS) technique described by Hedley etal. in IEEE Trans. Med. Imag., Volume 10, 1991, pages 548 to 553 andextended by Gmitro et al. in “Information Processing in Medical Imaging,ed. Bizais, Barillot, Paola, Iles de Berder, Kluwer, 1995. This has beenproposed as a means of correcting for motion during diffusion weightedimaging of the brain.

In the implementation used by Gmitro, the constraint in the image domainis provided by an image acquired without any diffusion weightinggradient (the reference image). Correction is performed in hybrid space(x vs. K_(y)), under the constraint that the magnitude and higher-orderphase terms of each line are correct.

A binary mask is generated from this reference image (spatial domain),in which pixels inside the head are given the value 1, and pixelsoutside the head the value 0. This constraint is applied to thecorrupted images by multiplying the image by the mask function. Theresulting product image has a black background, and consequently lacksthe ghosting artefacts outside the head that are caused by motion. Thismasked image is then transformed into hybrid space, and is compared,line by line, with the hybrid space of the original (unmasked) image.The zero and first order phase terms that, when applied to the originalhybrid space lines, produce the closest match (in a root mean squaresense) to the lines of the hybrid space from the masked image are thenfound using a search procedure. These phase terms are then treated asestimates of the motion correction parameters, and the algorithmreiterates.

In general, a mask image as required by this algorithm, may not beavailable, or may be out of registration with the images that needcorrection, and thus inaccurate.

The use of averaging or correlation to smooth out motion artefacts hasbeen described in U.S. Pat. Nos. 5,363,044, 5,233,302, 5,124,649, and4,966,149. The correction of image data using motion information derivedfrom the image data or from extra motion detection or position trackingsequences is described in U.S. Pat. Nos. 5,427,101, 5,382,902, 5,254,948and 5,251,128 and in Japanese Pat No 05080248. U.S. Pat. No. 5,311,132describes correcting demodulation frequency in an MRI imaging processusing a focusing criterion to determine image blur.

The technique described in U.S. Pat. No. 5,363,044 requires theacquisition of two interleaved data sets. These data sets are combinedinto a single image. The image is divided into sub-sections and a“gradient energy” term is summed over all the pixels of each sub-sectionof this image. The gradient energy term is defined as the sum over allthe pixels of the total of the squares of four partial derivatives ofthe real and imaginary parts of a complex image, with the partialderivatives being the difference between an image and the image shiftedby one pixel. A phase term used in the data combination is varied untilthe gradient energy sum in each sub-section is minimised. This techniquedetermines a phase factor which locally removes ghosts in an image. Thephase factor depends on the interleave conditions and the harmonicnumber of the motion. It assumes that the motion is periodic orquasi-periodic with respect to the acquisition order. It is likely thatto be effective, this technique requires the image to be sub-dividedinto regions. where the ghosts do not overlap either each other or themain body of the image and therefore its applicability is limited.

U.S. Pat. No. 5,311,132 describes a technique for correcting magneticresonance images by removing blur introduced as a result of magneticfield inhomogeneities and variations in magnetic susceptibility of anobject being imaged. It does not correct for artefacts due to motion.The technique described therein involves demodulating acquired data atdifferent frequencies to overcome such blurring. Using a local focuscriterion, the demodulation frequencies which optimise the focus of eachregion of an image are determined. Compensating for magnetic fieldinhomogeneities is significantly less complex than correcting for motionartefacts. Only one parameter needs to be searched for at each imageregion. Correction for motion artefacts is more complex because motionaffects the focus in all regions and in general there are six degrees offreedom rather than one for field inhomogeneities. The technique isapplicable only to non-Fourier transform reconstruction techniques.

SUMMARY OF THE INVENTION

It is an aim of the present invention to provide an alternative methodfor processing signals generated by a magnetic resonance imaging machineto generate an image such that image defects arising from patient motionduring the acquisition of the signals are corrected for in a mannerwhich does not require a modified signal acquisition technique.

The present invention provides a method for generating an artefactreduced physiological image of an object comprising the steps of:

(i) acquiring a data set from a plurality of received signals generatedin an object measuring process;

(ii) manipulating the data set in order to generate an image therefrom,wherein the data set manipulation includes the reduction of imageartefacts resulting from object motion during said object measuringprocess,

characterised in that the data set manipulation step includes the stagesof:

(a) calculating a focus criterion for an initial artefact containingimage;

(b) generating an initial model for the motion of the object during theobject measuring process;

(c) manipulating the data set in order to compensate for the effects ofthe model and recalculating the focus criterion for the manipulated dataset; and

(d) iteratively varying the model and repeating stage (c) in order toobtain a final object motion model with which the focus criterion isoptimised;

wherein the artefact reduced image is then generated from the data setafter the data set has been manipulated in order to offset the effectsof the final object motion model.

The invention provides the benefit of providing a method whereby theimage quality of certain images which are distorted by motion of theobject during data acquisition may be improved.

The method of the invention is particularly applicable for manipulatingimages obtained using magnetic resonance imaging (MRI) techniques.Preferably the focus criterion is an image entropy criterion. The termimage entropy is taken here to mean the degree of disorder in an image.

Patient motion in MRI scans may be in directions known as the readoutdirection and the phase encode direction. The data set may bemanipulated to reduce the effect of patient motion in these directions.In addition, rotation motions may be corrected.

An image entropy criterion may be taken over the whole of an image or itmay be determined for specific regions of an image. For example in anMRI image of a head, the image entropy may be calculated for only forthose regions of the whole image in which the patients head is imaged.

Several techniques for varying the model of patient motion may be used.The patient motion model may be varied in a piecewise linear orpiecewise constant manner. In more sophisticated techniques segments ofthe patient motion curve are varied with the size of the segments beingreduced as the model is optimised. Alternatively a limited set of nodesof the patient motion curve may be initially varied with more nodes thenbeing varied to obtain a finer scale patient motion curve. It may bedesired to initially obtain an image over a limited sub-set of k-spacelines, with the focus criterion for this sub-image being optimisedbefore the sub-set is expanded. These techniques might be employed toreduce the time taken to obtain a corrected image.

In a further aspect, the invention provides a method of improving theimage quality of magnetic resonance imaging images by the reduction ofmotion induced artefacts comprising the steps of:

(a) calculating a focus criterion for an initial artefact containingimage;

(b) generating a model of a possible motion sequence as a cause of theartefact;

(c) manipulating the data set in order to compensate for the effects ofthe model and recalculating the focus criterion for the manipulated dataset; and

(d) iteratively varying the model in order to obtain a final model withwhich the focus criterion is optimised;

wherein the artefact reduced image is then generated from the data setafter the data set has been manipulated in order to offset the effectsof the final model.

Compared with prior art techniques for MRI image improvement where aspecialised image acquisition strategy has to be adopted, this aspect ofthe invention provides the advantage that a conventional imaging processmay be used which is then followed by a post-acquisition processingstage.

In another aspect, the invention provides a magnetic resonance imagingscanner arranged to generate an artefact reduced image of a patient byperforming the steps of:

(i) acquiring a data set from a plurality of received signals generatedin a patient measuring process;

(ii) manipulating the data set in order to generate an image therefrom,wherein the data set manipulation includes the reduction of imageartefacts resulting from patient motion during data acquisition,

characterised in that the data set manipulation step includes the stagesof:

(a) calculating a focus criterion for an initial motion artefactcontaining image;

(b) generating a patient motion model as a possible cause of theartefact;

(c) manipulating the data set in order to compensate for the effects ofthe model and recalculating the focus criterion for the manipulated dataset; and

(d) iteratively varying the model in order to obtain a final model withwhich the focus criterion is optimised;

wherein the artefact reduced image is then generated from the data setafter the data set has been manipulated in order to offset the effectsof the final model.

Subsequent to the filing of the priority application, a scientific paperby R. A. Zoroofi, Y. Sato, S. Tamura and H. Naito was published in IEEETransactions on Medical Imaging, Vol. 15, No. 6, December 1996, pages768-784, relating to MRI artefact cancellation due to rigid motion inthe imaging plane. The paper concerns blurring and ghosting due tomotion during the image acquisition. The technique described thereininvolves using a minimum energy method to estimate unknown motionparameters. The minimum energy method requires a knowledge of theboundary of a region of interest whereas an entropy focus criterion doesnot require knowledge of any boundary position. The method of Zoroofi etal. is to select a group of k-space lines over which they believe therehas been no motion, based on tests which are performed algorithmically,and an image is formed from these lines, thresholded and used as theboundary of the region of interest. In order for the method to beeffective it is necessary to be able to identify a group of k-spacelines over which there is little motion. The method of the presentinvention does not have this requirement and can compensate forcontinuous motion.

BRIEF DESCRIPTION OF THE DRAWINGS

In order that the invention may be more fully understood, examplesthereof will now be described, by way of example only, with reference tothe accompanying drawings in which:

FIG. 1 shows a schematic diagram of a magnetic resonance imaging system;

FIG. 2 is a functional diagram of the operation of the FIG. 1 system;

FIG. 3 is a flow diagram of an image correction routine performed by theFIG. 1 system;

FIG. 4 is a example of a hypothetical patient motion curve;

FIG. 5 shows a representation of the perturbation of a patient motioncurve;

FIG. 6 shows MRI images before and after motion correction for motion inthe readout direction;

FIG. 7 shows MRI images before and after motion correction for motion inthe phase encode direction;

FIGS. 8 and 9 are pseudo-code algorithms; and

FIGS. 10 to 12 are schematic illustrations of further image correctiontechniques;

DETAILED DISCUSSION OF PREFERRED EMBODIMENTS

Referring to FIG. 1 there is shown a schematic diagram of a magneticresonance imaging system 10. The system 10 incorporates a magneticresonance imaging scanner 12 of conventional type. The scanner 12 has asuperconducting main magnet 20 which generates a magnetic fieldsufficiently strong to cause a net alignment along the field directionof atomic nuclei within a patient. The scanner 12 also includes shimcoils 22 in order to correct for undesired inhomogeneities in themagnetic field of the main magnet 20. The magnetic field produced by theshim coils 22 is controlled by a shim coil power supply unit 24.

The resonance frequency of particular atomic nuclei is characteristic ofthe nucleus and the strength of the applied magnetic field. In order toprovide spatial information, a magnetic field gradient is generated bygradient coils such as coils 26. Gradient coils are often arranged togenerate gradient fields in three orthogonal directions. The magneticfields generated by the gradient coils are controlled by a gradient coilpower supply unit 28. In order to generate a signal from the atomicnuclei of the patient, a radio-frequency magnetic pulse is generated bya transmit coil 30. This pulse “flips” the angle of the nuclear spinswithin a certain patient slice of volume. These “excited” spins, ormagnetisations, then induce a current in a receive coil which may be thesame coil as the transmit coil 30. The coil 30 is connected to atransmit unit 32 and a receive unit 34, each of which also receivessignals from a frequency source 36.

The system 10 includes a controlling computer 38 which controls theoperation of the components of the system 10. The computer 38 controlsthe gradient coil power supply unit 28 in the form of gradient timing,magnetic field strength and orientation control. In addition, thecomputer receives signals from the receive unit 34 together withtransmitter timings.

In order to form an image of the organs of a patient, the patient isinserted into the system 10 and a series of measurements are taken withdifferent combinations of static and/or varying gradient fields. Thesignals from the tissue of the patient depend on the tissue'sproperties, the magnetic field gradient strengths, gradient orientationsand timings with respect to the applied radio frequency pulses. Thevarying gradients code the received signal's phase, frequency andintensity. The received signals as a function of time form an orderedset which is stored in memory in the computer 38 for subsequentprocessing.

In the subsequent signal processing stage, a Fourier transform isperformed on the ordered set of received signals, with the modulus ofthe transform being used to assign the signals to a grey scale in orderto form an image. The set of received signals is often referred to asexisting in “k-space”.

In conventional MRI, if the patient moves during the acquisition ofdata, the received signal is affected and part of the k-space signal iscorrupted. Because of the way the image is reconstructed, this motionaffects the whole image, causing blurring and/or ghosting artefacts inthe final image.

Referring now to FIG. 2 there is shown a functional block diagram of theoperation of the system 10. The computer 38 controls and receivesinformation from the scanner 12 and uses this information to generate animage on a display 50. This image is an initial reconstructed image. Ifan operator of the system 10 considers that the initial image iscorrupted due to patient motion, a further signal processing routine isselected. Data which is stored in the computer 38 is then furtherprocessed in order to reduce the effects of the patient's movement. Thefollowing is a description of the principles behind this further signalprocessing routine.

The magnet, or gradient frame is denoted by capital letters, the bodyframe by lower case. Consider a point P in the gradient frame. Itsposition in the body frame, is {right arrow over (r)}₀ before movementand {right arrow over (r)} after movement. In the gradient frame, theposition of P remains at {right arrow over (R)}, irrespective of thebody's motion. At some pre-defined time, for example at the middle ofthe acquisition, when k_(y)=0, the body and gradient frame origins aretaken to coincide and {right arrow over (r)}₀={right arrow over (R)}.The patient motion with respect to this time and the gradient frame isdescribed by a vector {right arrow over (M)}. Before and after motion,the magnetisations in the gradient frame ρ_(G) are related to the bodyframe magnetisations ρ_(b) (the desired image) by: $\begin{matrix}{{\rho_{G}\left( \overset{\rightarrow}{R} \right)} = \quad {\rho_{b}\left( {{\overset{\rightarrow}{r}}_{0} = \overset{\rightarrow}{R}} \right)}} & {\quad \left( {{before}\quad {motion}} \right)} \\{= \quad {\rho_{b}\left( {\overset{\rightarrow}{R} - \overset{\rightarrow}{M}} \right)}} & {\quad \left( {{after}\quad {motion}} \right)} \\{= \quad {{\rho_{b}\left( \overset{\rightarrow}{r} \right)}.}} & \quad\end{matrix}$

The signal S (ignoring constants, relaxation etc.) is measured in thegradient frame

S({right arrow over (K)})=∫ρ_(G)({right arrow over (R)})exp[i{rightarrow over (K)}.{right arrow over (R)}]d{right arrow over (R)}

where K is the k-space co-ordinate in the gradient frame. So aftermotion

 S(K)=∫ρ_(b)({right arrow over (R)}−{right arrow over (M)})exp[i{rightarrow over (K)}.{right arrow over (R)}]d{right arrow over (R)}

ie

S(K)=∫ρ_(b)({right arrow over (R)}−{right arrow over (M)})exp[i{rightarrow over (K)}.({right arrow over (R)}−{right arrow over (M)})]exp[i{right arrow over (K)}.{right arrow over (M)}]d {right arrow over (R)}

It is possible to swap d{right arrow over (R)} for d{right arrow over(r)} because {right arrow over (r)}={right arrow over (R)}−{right arrowover (M)} and does not depend on {right arrow over (R)} (for a rigidbody, all parts of the body move together independent of position in thegradient frame). Hence

S({right arrow over (K)})=∫ρ_(b)({right arrow over (r)})exp[i{rightarrow over (K)}.({right arrow over (r)})]exp[i{right arrow over(K)}.{right arrow over (M)}]d{right arrow over (r)}.

The motion vector {right arrow over (M)} is composed of a translation{right arrow over (M)}_(T) and a rotation {right arrow over (M)}_(R).Any rigid-body motion can be described as a translation followed by arotation. Here the rotation axis is taken to be through the origin ofthe body frame and rotation is described by the small angleapproximation {right arrow over (M)}_(R)={right arrow over (θ)}×{rightarrow over (r)}.${\overset{\rightarrow}{M}}_{R} = {{\overset{\rightarrow}{\theta} \times \overset{\rightarrow}{r}} = {{\det {\begin{matrix}\overset{\rightarrow}{i} & \overset{\rightarrow}{j} & \overset{\rightarrow}{k} \\0 & 0 & \theta_{z} \\r_{x} & r_{y} & r_{z}\end{matrix}}} = {{\overset{\rightarrow}{i}\left( {{- \theta_{z}}r_{y}} \right)} + {\overset{\rightarrow}{j}\left( {\theta_{z}r_{x}} \right)}}}}$

for θ_(x)=θ_(y)=0. This cross product is the change in co-ordinates notthe new one.

The translation may be brought out of the integration as it has no{right arrow over (r)} dependence (and in the absence of rotation wehave the Fourier shift theorem).

S({right arrow over (K)})=exp[i{right arrow over (K)}.{right arrow over(M)} _(T)]∫ρ_(b)({right arrow over (r)})exp[(i{right arrow over(K)}.{right arrow over (r)}]exp[i{right arrow over (K)}.{right arrowover (M)} _(R) ]{right arrow over (r)}=

exp[i{right arrow over (K)}.{right arrow over (M)} _(T)]∫ρ_(b)({rightarrow over (r)})exp[i{right arrow over (K)}.{right arrow over(r)}]exp[i{right arrow over (K)}.({right arrow over (θ)}×{right arrowover (r)})]d{right arrow over (r)}

Cyclic permutation of the scalar triple product gives {right arrow over(r)}.({right arrow over (K)}×{right arrow over (θ)}) and hence

 S({right arrow over (K)})=exp[ i{right arrow over (K)}.{right arrowover (M)} _(T)]∫ρ_(b)({right arrow over (r)})exp[i{right arrow over(r)}.( {right arrow over (K)}+d{right arrow over (K)})]d{right arrowover (r)}

With rotation only about the z axis, d{right arrow over (K)}={rightarrow over (i)}({right arrow over (K)}_(y){right arrow over(θ)}_(z))−{right arrow over (j)}(K_(x)θ_(z)).

This is the Fourier rotation theorem with a small angle approximation.For image correction, translations can be corrected by applying thephase-shift factor exp[−i{right arrow over (K)}.{right arrow over(M)}_(T)] to the measured signal. To compensate for rotation, eitherk-space or the image should be rotated. For diffusion weighted images,the effect of rotation is to produce a d{right arrow over (K)} thatdepends only on the diffusion gradients, not the imaging gradientsrepresented by K_(x) and K_(y). This makes the compensation problemeasier as rotation causes a constant shift of k-space and can becorrected by a phase-shift in the image domain.

During the measurement process, N_(y) readout lines are obtained. Ifthere is patient motion during the measurement process then there is aresulting displacement of the patient at each readout time. The N_(y)readout lines are indexed by the superscript j where j=0,1,2 . . . ,N_(y)−1. Considering translational motion only, a displacement vector{right arrow over (M)}_(T) ^(j) at each time is composed of adisplacement m_(x) ^(j) in the frequency encoding direction (x) anddisplacement m_(y) ^(j) in the phase encoding direction (y) such that{right arrow over (M)}_(T)={right arrow over (i)}m_(x)+{right arrow over(j)}m_(y), dropping the superscripts for clarity. As given above, imagecompensation is performed by multiplying the k-space data by the phaseshift factor exp[−i{right arrow over (K)}.{right arrow over (M)}_(T)^(j)]. The compensation for translation is a linear process, as is theFourier transform between image and k-space. Consequently, correctionfor two dimensional translations can be separated into one dimensionalreadout and phase encode compensations.

For one dimensional readout compensations the phase term used to correcteach point u in the readout line (in k-space) is exp[−iK_(x).m_(x) ^(j)]where$K_{x} = {\frac{2\pi}{{FOV}_{x}}\left( {u - \frac{N_{x}}{2}} \right)}$

and FOV_(x) denotes the field of view in this direction and is the sizeof the imaged area, which typically is 250 mm×250 mm. This indexing usesthe convention that K_(x) increases from negative to positive as uincreases in steps of 1 from 0 to N_(x)−1, there being N_(x) data valuesin the readout direction.

For the one dimensional phase encode compensation, at each readout linethe correction for the displacement m_(y) ^(j) in the phase encodedirection is exp[−iK_(y).m_(y) ^(j)] where$K_{y} = {\frac{2\pi}{{FOV}_{y}}{\left( {j - \frac{N_{y}}{2}} \right).}}$

For correction due to motion in the phase encoding direction, everypoint in one readout line is multiplied by the same complex factor,whereas for motion in the readout direction, the phase factor variesthrough K_(x) at each point along a line.

In the image correction signal processing routine, the motion of thepatient is parametised into a finite number of discrete movements at anumber of sample points referred to as nodes. Whilst a plot of thepatient's motion against time is more realistically a smooth curve overthe imaging interval T, it is assumed that this motion may be wellapproximated by a piecewise constant, a piecewise linear or a cubicspline curve defined on the nodes. There may be one node for eachreadout line, alternatively one node could represent a group of severalconsecutive readout lines. If k-space is not traversed sequentially, theordering of nodes should represent the temporal order of theacquisition.

In the motion correction routine, in the absence of any prior knowledge,it is initially assumed that the patient was in fact stationary duringthe data acquisition. Perturbations are then introduced into thepatient's position at the nodes and the effect that these perturbationshave on the resulting image quantified. If the image is improved by sucha perturbation, as reflected by an improvement to a focus criterion,then this is indicative that the perturbation introduced in some wayreflects the motion of the patient. By iterating round a loop in whichthe perturbations are varied and analysing the effect this has on thefinal image, an approximation to the actual motion of the patient duringthe course of the data acquisition is obtained, and this information isused to correct the image in a manner described above, where a linearmotion is corrected for by multiplying each k-space line by a respectivephase factor.

The preferred measure of the image “quality” is an entropy focuscriterion. The term entropy is taken here to mean the degree of disorderin an image. A conventional measure of image entropy, E, is given by theequation:${E = {\sum\limits_{x,y}\quad {\frac{I_{x,y}}{I^{\max}}\log \frac{I_{x,y}}{I^{\max}}}}},$

where I_(x,y) is the image intensity at a point x,y and I^(max) is themaximum intensity in the image. Here intensity is defined as the modulusof the complex value at a particular point. The denominator I^(max) mustnot be less than any value of I_(x,y) at any point x,y upon anyiteration. Phase compensation in MRI involves the addition of manycomplex numbers, with the final image being the detected magnitude ofsuch sums, and this depends rather sensitively on the exact phasecompensation applied. As a result the image maximum is variable,depending on the exact estimate of patient motion applied to calculatethe image quality. For example, it may happen that the image maximum fora fully corrected image is much lower, perhaps by a factor two, than inthe initial uncorrected image used as a starting point. A low maximumcauses the above measure of global entropy to be lowered, which wouldnot be desirable. Therefore a modified entropy focus criterion:${E = {\sum\limits_{x,y}\quad {\frac{I_{x,y}}{I^{maj}}\log \frac{I_{x,y}}{I^{maj}}}}},$

is more suitable for motion correction, where I^(maj) is some numberwhich is larger than any perturbed image maximum likely to beencountered. For an uncorrected image J₀, if its modulus I₀ is taken,this will have some maximum I₀ ^(max). If this value is doubled, asuitable value for I^(maj) is obtained for all images J derived from theimage J₀.

Referring now to FIG. 3, there is shown a flow chart 100 of stages in animage processing routine for reducing the effect of patient motion inone dimension. The image processing is carried out on previouslyacquired data, as indicated by a box 102, obtained as sequential k-spacelines, having an order starting at low negative frequencies andprogressing through DC at a time t⁰ to high positive frequencies. Ifnecessary, the data in k-space may be reordered as indicated by a box104. FIG. 3 shows the stages in a routine which performs the functionsof:

i) generating an initial patient motion spline parameterisation;

ii) implementing a set of spline perturbations on the current bestestimate of the motion spline;

iii) using the perturbed motion spline to obtain a compensated image;

iv) calculating a focus criterion for the compensated image;

v) comparing the focus criterion of the compensated image with that ofthe image prior to the perturbation.

In the flow chart 100, acquired data is passed to a patient motion curveinitialisation stage 106. FIG. 4 shows a graph 108 showing a plot 110 ofthe movement of a patient during the image acquisition. The y-axis ofthe graph 108 is patient displacement and the x-axis represents time.Whilst the patient's motion is continuous, image acquisition shots aretaken at discrete moments in time, as indicated by dots 112. Thepatient's motion is represented by a set of nodes, N, which may be lessthan the number of image acquisition shots. The initialisation stage 106sets the initial patient motion curve to one in which the patient isstationary during the whole measurement period and hence the patientmotion curve would lie along the x-axis of the graph 108.

In a node selection stage 114, a node n which is a member of the set ofnodes N is selected, nεN. A perturbation to the node n is selected andin a perturbation stage 116, the node n is moved, or perturbed, by adisplacement parallel to the y-axis of graph 108 by the selectedperturbation. Initially the perturbation will have a step size of 2 mm.If the motion curve of the patient is considered to be an N-dimensionvector {right arrow over (P)}, then after the perturbation the vector is{right arrow over (P)}+{right arrow over (δ)} where {right arrow over(δ)} is from a perturbation set Δ_(n), {right arrow over (δ)}εΔ_(n), andin general is dependent on n.

FIG. 5 shows a schematic representation of the effect on a patientmotion curve of the perturbation stage 116. FIG. 5 shows arepresentation of the perturbation stage on an initialised patientmotion curve and the effect of the stage on an intermediate patientmotion curve. In FIG. 5, prior to the perturbation stage, the patientmotion curve is represented by a solid line 120. The perturbation stagehas the effect of adjusting the motion curve at a discrete time point,as indicated by a dashed line 121. The intermediate patient motion curvebefore the perturbation stage is represented by a solid line 122 withthe perturbation stage having the effect of adjusting the intermediatemotion curve at a discrete time point as indicated by a dashed line 123.

Referring again to FIG. 3, after the perturbation stage 116, the imageprocessing routine applies the perturbation of this stage in a motioncompensated data formation stage 124. In the stage 124, the k-space datais manipulated by the application of a phase correction, as describedpreviously, with the phase shift to be applied being dependent on theperturbation to the patient motion curve. Once the k-space data has beenmanipulated, a Fourier transform stage 126 performs a 2-D fast Fouriertransform on the manipulated data. In this stage the modulus of thecomplex Fourier transformJ_({right arrow over (P)}+{right arrow over (δ)}) is also obtained,generating a modified imageI_({right arrow over (P)}+{right arrow over (δ)}). The focus criterion,F(I_({right arrow over (P)}+{right arrow over (δ)})) δεΔ, for thismodified image is then calculated using the modified entropy focuscriterion as given above. In a comparison stage 130, the focus criterionF(I_({right arrow over (P)}+{right arrow over (δ)})) δεΔ is comparedwith the focus criterion of the unperturbed image,F(I_({right arrow over (P)})).

IfF(I_({right arrow over (P)}+{right arrow over (δ)}))<F(I_({right arrow over (P)}))δεΔ then the perturbation δ has a beneficial effect on the overall imageand so the motion vector {right arrow over (P)} is updated so that{right arrow over (P)}={right arrow over (P)}+{right arrow over (δ)} bya motion vector update stage 132, otherwise {right arrow over (P)}remains unchanged. The routine cycles round in a loop updating {rightarrow over (P)} for the whole set of nodes until the decrease in entropyon traversing the whole node set is less than a set value, typically0.2. The step size of the perturbation is then halved to 1 mm to give afiner splice approximation and the process repeated, followed by afurther halving of the perturbation magnitude to 0.5 mm. Once theiterative process has been completed for this step size, the finalpatient motion vector is used to modify the data set from which a fullycorrected image is obtained. The fully corrected image is then storedand displayed on a screen for analysis by a medical practitioner.

Referring now to FIG. 6 there are shown saggital MRI images before andafter motion correction. The readout direction is head to foot and thephase encoding direction is anterio:posterior. FIG. 6a shows an imagewhich has been corrupted by periodic patient motion along the readoutdirection. FIG. 6b shows an image obtained from the FIG. 6a data afterthe motion correction routine has modified the data.

FIG. 7 shows saggital images before and after correction for motion inthe phase-encode direction. FIG. 7a shows an image corrupted by a stepmotion during image acquisition. FIG. 7b shows an image after correctionof motion using knowledge of the time and measure of the displacementand FIG. 7c shows an image after the correction routine has compensatedfor the motion. The FIG. 7c image is not as good as the manually focusedimage of FIG. 7b, with the manually focused image having a lower entropythan that of FIG. 7c, although the entropy of the FIG. 7c image is lowerthan that of the FIG. 7a image. It is possible that the image correctionroutine has found a solution which lies in a local minimum of theentropy function for the given set of perturbations.

Referring to FIG. 8 there is shown a pseudo-code listing of an algorithmto perform the stages of the FIG. 3 routine for data acquired in a dataacquisition process in which the k-space lines have been obtainedsequentially. Prior to the execution of the algorithm, the nodes in thepatient motion spline, called m(t) are ordered {1, 2, . . . , N} where Nis the number of nodes of the set of nodes N and is usually an oddnumber. The middle node is denoted midNode=(N−1)/2. A node orderingpermutation π(.) is defined as π(1)=midNode, π(i)=midNode+i/2 if i iseven, π(i)=midNode—(i−1)/2 if i is odd for i ε (1,2, . . . , N). Thisnode ordering corresponds to working out radially from the mid k-spacenode which for most common acquisitions is at DC or zero frequency.Perturbing nodes near the mid-node tends to have the biggest effect onreducing entropy. An algorithm termination threshold, ε, is set at avalue of for example 0.2. In general the centric node ordering is notmandatory and would not be used for all acquisition modes.

Whilst the above description is a description of a routine for reducingthe effects of motion of a patient in one direction, it may be possibleto extend the principles of the routine to correct for movement in twodimensions, including rotational movements. Instead of a 1-D motionfunction m(t) to identify, in the 2-D case there are threem(t)=m_(x)(t), m_(y)(t), m_(r)(t) where the motion curve m_(r)(t) is therotation. Using the same number and spacing of nodes as before, thediscretised parameter space, which was described above, is now 3 timesas large as for the 1-D case. The node ordering will again be centricfrom the DC line in the middle of k-space. At each node in a givensweep, the question arises as to what order to thread through thefunctions m_(x)(t), m_(y)(t), m_(r)(t)=m₁(t), m₂(t), m₃(t) performing a‘node jiggle’ on each in turn at the given node. This could bedetermined by experimentation—it is unlikely to be critical, but as ageneral rule it is best to put that parameter likely to give the biggestdescent (entropy reduction) first, as this will give the fastest andmost robust convergence. This was also the reason the centric nodeordering described above was chosen. For the 2-D case, the majormodification to the above pseudo-code is thus to insert a further doloop immediately nested inside the do 200 loop above, as shown in FIG.9.

An alternative possibility for 2-D motion correction would be to reducethe dimensionality of the problem by rendering the three motionparameters at least partly separable. Taking the middle of the image asthe nominal rotation centre O, then rotation has little defocusingeffect in the neighbourhood of the centre. Thus a focus criteriondefined on a limited rectangular region circumscribed about O could inprinciple be employed to estimate the translational motion in theneighbourhood of O only. Once the translation has been estimated andcorrected for, a similar procedure could be used to estimate therotation separately. This time the focus criterion would be evaluated onselected outer regions where any rotational blurring would be mostpronounced.

Rather than using an entropy expression which is defined globally over awhole image, entropy expressions may be defined in a more localisedmanner. Instead of calculating the entropy over the whole image, it maybe calculated “linewise”, that is within an x-resolution cell ory-resolution cell respectively. This would be useful if it were knownthat the patient motion was confined to the y- or x-directionsrespectively. This would not therefore be appropriate for the general2-D patient motion case. The modified linewise entropy expression may begiven as:$E = {\sum\limits_{x}\quad {\sum\limits_{y}\quad {\frac{I_{x,y}}{I_{x}^{maj}}\log {\frac{I_{x,y}}{I_{x}^{maj}}.}}}}$

In this expression it is better to include in the summation overx-resolution cells only those cells whose intersection with the initialimage contains substantial energy over a certain threshold, otherwiseoff-image lines containing only defocus energy could be problematical.

It may be beneficial to calculate the entropy not on the whole image butrather on selected sub-regions of the image, for example rectangularregions. The sum over all or selected sub-region entropies could then beformed to give a global image focus criterion. Care would need to beexercised in the implementation of this concept since image “mass” orintensity can migrate from one sub-region to another due to blurring,and thus one sub-region is not independent from its neighbours. Theeffect of this migration is less if the sub-regions are taken to bequite large rectangles, occupying at least a few percent of the totalimage area.

Ghosting in images tends to occur along the phase encode direction. Theimage resolution may be reduced in the readout direction withoutaffecting the node placements or the resolution of ghosts in the phaseencode direction. The focusing can be applied on the reduced resolutionimage to obtain the motion spline. This is then applied retrospectivelyto a full resolution original image to obtain a full resolution motioncorrected image. This technique has the benefit of reducing computationtime since the reduced resolution image contains fewer points.

An alternative strategy for node perturbation is to allow theperturbation of each node to take any arbitrary value to minimise theentropy, rather than limiting the perturbation to a few integralmultiples of 2 mm, for example.

Entropy-based focus criteria are not the only possible choices for afocus criterion. It would be preferable to have a focus criterion whichcould be computed directly in k-space or even hybrid-space, since thiswould result in a considerable saving in computing time. The focuscriterion is a function which is optimised when the patient motion iscorrectly estimated, either being maximised or minimised at this point.

As an example of an alternative type of focus criterion, an imagevariability focus criterion might be employed, where the focus measuredetects the presence of a large number of oscillatory terms in adefocused image. The measure of image variability is a standardstatistical one: (image standard deviation)/(mean intensity) or(variance)/(mean squared intensity). The latter may be expressed asfollows:${{F_{IV}(J)} = \frac{\sum\limits_{x,y}\quad \left( {I_{x,y}^{2} - \overset{\rightarrow}{I}} \right)^{2}}{{\overset{\rightarrow}{I}}^{2}}},$

where {right arrow over (I)} is the image mean squared ΣI_(x,y) ².I_(x,y) ² is the magnitude squared (or intensity) of the x,y^(th) pixelin the modulus image I. If at focus the image is characterised byminimum variability, then this focus criterion will be minimised atfocus. The normalisation by the mean is designed to compensate forvariations in the absolute image level. This is not strictly necessarysince by Parseval's theorem of

Fourier analysis, {right arrow over (I)} for a phase-corrected imagewill be equal to {right arrow over (I)} before any compensation. Thusthis last equation may be simplified to:${F_{IV}(J)} = {\sum\limits_{x,y}\quad {I_{x,y}^{4}.}}$

The operation of the motion correction routine may be modified in avariety of ways. For example, the node set N and the sequence in whichnodes are traversed may be altered and may be altered dynamically as thealgorithm proceeds. The perturbation set Δ_(n) may be modified, as maythe type of spline approximant. The optimisation algorithm may beadapted, as may the exact manner in which the image compensation isperformed. For example rather than perturbing a single node, it may bebeneficial to perturb two or more nodes together, with the algorithmvarying the amplitude of each node perturbation.

In fast scanning MRI image acquisition, such as turbo spin-echo imaging,k-space is rapidly traversed using only a limited number of equallyspaced lines in the phase-encode direction. It is then traversed again,moving along by one line until all of k-space has been acquired. Eachsuccessive sub-image is formed from slightly different k-space lines,and only one will contain the DC line. As the patient moves, successivesub-images will be translated, and, in the 2-D generalisation rotated,relative to each other reflecting this movement. It may be possible todetermine such movement approximately by maximising the normalisedcross-correlation of adjacent sub-images. This method may be a usefulstarting point for a further routine using the previous focus criterionapproach. It may also be possible to use the focus criterion approach,searching through a number of possible translational or rotationalshifts to optimise a focus criterion on the composite sub-image. It maybe possible to merge the k-space lines of sub-image pairs which areinterlaced in k-space to form a composite sub-image. This may then berepeated on the higher resolution sub-images formed, with the relativetranslation and rotation of these sub-images being estimated via across-correlation. This merging could continue until a final image isformed, encompassing the whole of k-space.

A benefit of the above image processing techniques is that the dataacquisition stage is unaffected and therefore the technique may beapplied retrospectively to previously acquired images. Indeed the imageprocessing may be carried out on a computer which is independent of theMRI scanner.

Rotational motion part way through a scan causes a rotation of thecorresponding section of k-space. If the image is to be reconstructedfrom k-space using an FFT, the rotated data should be interpolated or“re-gridded” to lie on a regular Cartesian grid. Direct integration withan irregular grid can be used instead of the FFT but this is slower. Analternative option is to perform rotation compensation in the imagedomain but this requires the original data to have a global phasecorrection applied to make the k-space Hermitian. Furthermore, thecorrection scheme requires the rotation to be discrete with periods ofconstant position between motion. Re-gridding is the preferred option.

A convolution kernel is chosen that is a balance between computationalspeed, accuracy and the unwanted introduction of extra image artefacts.A separable windowed sinc kernel h is taken, where

h(k)=[1+cos(πk/k _(r))]sinc(4k/k _(r)),

where k is the number of k-space points along the frequency encode orphase encode direction from the k-space position being interpolated tothe data point and k_(r), the kernel “radius” is 4.78, following atechnique described by S. R. Marschener and R. J. Lobb in Proc.Visualization '94, R. D. Bergeron and A. E. Kaufman, Ed., IEEE ComputerSociety Press, 1994, pages 100-107. Because the focus criteriontechnique is sensitive to ghosts, data is oversampled to reducepost-aliasing effects. The first sidelobes of the kernel transform thatcause post-aliasing into the original image region are 31 dB below thecentral region, which itself varies by no more than 0.5% of its value atthe image centre.

The convolution process sums data points from a region of k-space aroundthe position to be interpolated, applying the kernel weighting to eachpoint. Points in the summation are ignored if they are beyond the edgeof k-space or in the oversampled regions or undersampled regions thatoccur on rotation. To account for the un-even sampling, the kernelvalues h(k) are summed and this sum is used to normalise theinterpolated point. If this sum is too small (taken to be less thanh(k=0)/5) the interpolated value is set to zero. For rotation angles ofless than 15°, zero filling has been shown to be acceptable.

Rotation of a complete image should not change the image entropy. If theentire k-space of a good image is rotated by up to 15°, the largestentropy change from the original image is approximately 0.5% of theoriginal image entropy, which compares with normal entropy changes dueto motion compensation of approximately 2%.

This rotation of k-space rotates the image about its centre. For mostimages, the true centre of rotation will not be the image centre. Tocounteract this, the rotated sections of k-space must have an unknowntranslational correction applied after the rotation correction.

The real and imaginary data are separately interpolated. To overcomeproblems if the image is centred in the field of view because of the FFTprocess resulting in very rapid sign changes of the real and imaginaryk-space data, the image is shifted half a field of view in x and y withwrap around allowed, prior to transforming to k-space.

Trials using simulated motion on a good quality image with a nod-typemotion of between 0° and 5° within ±10 nodes of the k-space centre showthat it is possible to determine rotation to the nearest degree bytrying each of the combinations of five rotation angles and twenty steptimes followed by a simple gradient descent algorithm for thetranslational correction. The process used to compensate for rotationalmotion is to repeat the following procedure for each possiblecombination of rotation angle and step time:

(i) zero pad image to four times original area

(ii) shift image by half the new FOV in x and y allowing wrap around

(iii) perform FFT to k-space

(iv) apply re-gridding to compensate for trial rotation angle and steptime

(v) perform iFFT to image domain

(vi) shift image by half the FOV to correctly recentre it

(vii) extract the image (i.e. discard the previously zero-padded regionso that the new image is now the same size as the original)

(viii) apply a gradient descent algorithm to find the translationalcorrection for unknown centre of rotation.

Further distinct algorithmic methods have been developed and found to besuccessful on a large number of images. These further methods could beconsidered to be evolutionary developments of the “node jiggle” approachdescribed above. The methods differ in the type of successiveapproximation to the true patient motion and the type of search scheme.Normally, groups of k-space lines are perturbed rather than one at atime. For solving a general two dimensional problem of patienttranslations and rotations, the search space is potentially very largeand so it is beneficial to introduce short cuts which reduce the timetaken for a search and minimise the risk of the search becoming stuck ina local entropy minimum.

The first of these further methods is termed “Focus via SuccessiveApproximation” or FOSA for short, which is illustrated schematically inFIG. 10. In this method, the search routine starts with a large scalesegment length 200 of length 16 or 32 k-space lines. The position ofthis segment length 200 is varied to the right of the origin 202 and ateach position its height, or magnitude of the displacement, is varied inthe search for a large focus criterion descent. If such a large descentis found, the patient motion curve is updated. A similar search is thencarried out to the left of the origin 202 using a segment length 204.This combined process will yield the best piecewise constantapproximation to the patient motion at the crude scale used (16 or 32lines). The process is then repeated with the segment length reduced byhalf and again until a segment length of two lines is reached. In FIG.10, segment lengths 200 and 204 were found in the first phase, segmentlengths 206, 208, and 210 in a second phase and segment lengths 212 and214 in a third phase, approximating to an actual patient motion curve216. The additional segment lengths may effectively add to or subtractfrom displacements from earlier segment lengths, as illustrated bysegment lengths 212 and 214 respectively.

For the general two dimensional case, the search is in effect fivedimensional over rotation, phase encode and frequency encode amplitudesand segment position and length. This is potentially slow to achievecomputationally. It may therefore be preferable to do a search onrotation amplitude first, since this is the most likely motion to occurfor cranial scans. Pure phase encode or frequency encode translationswould seldom occur in isolation without being accompanied by rotation.Thus if a rotation is detected, a search on phase encode or frequencyencode translations can then be performed. Having determined phaseencode or frequency encode displacements, the search can be alternatedback to rotation to possibly improve the first estimate.

The second further method is termed “Advanced Node Jiggle” or ANJI forshort and is illustrated schematically in FIG. 11. As with the FOSAmethod, the ANJI method progresses from a crude large approximation tothe patient motion to a smaller scale approximation capable of followingthe finer details of movement. Initially four nodes 220 to 226 acrossk-space, separated by a segment length of thirty-two k-space lines(assuming a total of 128 k-space lines) are perturbed. Node perturbationis performed in pairs of nodes, with the nodes of each pair being onopposite sides of the origin 228. Thus nodes 220 and 222 form a pair, asdo nodes 224 and 226. For each pair, the optimum displacement isdetermined for reducing the focus criterion, with the patient motionbetween nodes being assumed to be piecewise linear. For each pair ofnodes, the amplitudes in the frequency encode and phase encodedirections are perturbed jointly. In a further stage, the separationbetween nodes is halved and the process repeated, proceeding outwardsfrom the centre of k-space radially. A given stage can be repeated ifdesired until there is no further focus criterion descent and the searchthen proceeds to the next stage until the segment length becomes four ortwo lines at which point the search terminates.

To extend the ANJI technique to fully two dimensional correction,encompassing patient rotation as well as translation, joint searchingwould be computationally prohibitive and therefore the search wouldalternate between nodes 220 and 222 and then nodes 224 and 226, andsimilar for further stages. As with the FOSA method, it may be possibleto restrict the main amplitude search to rotation, alternating withjoint phase encode and frequency encode searches at nodes where arotation is detected via the occurrence of a descent.

The third further method is termed “Increasing Phase Encode ResolutionSearch” or IPERS for short. The FOSA and ANJI methods developcompensated images and apply the focus criterion to these at the fullphase encode resolution, i.e. they use all the k-space lines indeveloping the image. It is however possible to develop an image over asubset of k-space lines, if the subset is symmetrical about the centreof k-space to avoid symmetry artefacts, and this “subimage” isself-contained in so far as it contains information on how to refocusthe image over its particular subset of lines, even if it may haverather low resolution in the phase encode dimension. There is noinformation in the higher frequency lines which could help determinepatient motion parameters associated with the lines in the subset. Thereis therefore no need to develop the image at the maximum resolution inorder to determine patient motion parameters near the centre of k-space.

The IPERS method is shown schematically in FIG. 12 and uses a piecewiseconstant approximation to the patient motion curve, with a segmentlength of typically four k-space lines. The method starts with twosegment lengths 230 and 232 on opposite sides of the origin 234. Segmentlength 230 is held at zero amplitude for the whole IPERS method, withpatient motions being determined relative to this reference. A search iscarried out on rotation amplitude and frequency encode and phase encodedisplacements on the segment length 232, looking for the largest descentin the focus criterion in the subimage formed on the current interval ofeight k-space lines only. On completion of this stage, the k-spaceinterval is extended by two segments 236 and 238 of length four oneither side of the origin 234, giving a k-space interval of sixteenlines on which to develop a subimage. The patient motion curve isinitialised at zero for the new segments and the amplitudes of these newsegments are then varied.

Rotation amplitude is first varied on the two segments then frequencyencode amplitudes then phase encode amplitudes. The patient motion curveis updated as usual when focus criterion descent occurs. The methodproceeds adding segments until the whole of k-space is covered, and thefull patient motion curve obtained.

Once patient motion is determined in a given region of k-space, it isdetermined once and for all. There is no need to revisit the same areaof k-space as the FOSA and ANJI methods do at varying scales.

As with the FOSA and ANJI methods, it is possible to alternate a givenamplitude search rather than to search jointly, thereby replacing a twodimensional search with effectively a one dimensional search and also toperform a frequency encode and a phase encode search only when a descenthas occurred in the rotation search. The entropy focus criteriondescribed above was calculated on the whole image. A Box Entropy FocusCriterion may be defined to make the focus criterion more sensitive tolower contrast images, for example brain tissue. The initial unfocussedimage I⁰ is divided into a number (typically sixty-four) equalrectangles, and the r maximum in each rectangle is determined. If themaximum is less than a set fraction, say 0.1, of the image globalmaximum I⁰ _(max), then the box is labelled “dark”. The “box entropy” ofa general image I is calculated as the sum over k of entropies BE^(k)calculated for each non-dark box. Thus if I_(i,j) ^(k) is the i,j^(th)pixel of the k^(th) box, and I^(0,k) _(max) is the maximum of the k^(th)box of the initial unfocussed image I⁰, the entropy of the k^(th) box isdefined as:

BE ^(k)=Σ(I _(i,j) ^(k) /I _(max) ^(0,k))log(I _(i,j) ^(k) /I _(max)^(0,k)).

In practice, a further scaling and thresholding of pixel values is doneto prevent low value pixels from overcontributing to the sum. It isdesirable to have a balance between high and low pixels in the image,with the contribution of the middle value pixels de-emphasised.

Whilst the invention has been described in relation to MRI imageimprovement, the hereinbefore described principles may be extendible toother imaging techniques where patient motion results in a degradationof the image quality, for example tomographic scanning using x-rays orultrasound scanning.

In addition to correction for patient motion, the techniques ofmodifying a data set in order to optimise a focus criterion may beapplicable to situations where timing errors or instrumental factorsresult in image degradation. For example, in EPI image acquisitiontiming errors result in Nyquist ghost images being formed. Suchartefacts may be corrected where there is a physical and mathematicalrelation between an image artefact and its physical cause.

What is claimed is:
 1. A method for generating a corrected image of anobject by: (i) employing scanning means (12) to scan an object andderive an image data set from signals received from the object, and (i)manipulating the data set to reduce unwanted artefact effects associatedwith the object and produce corrected image data to yield anartefact-reduced image, characterised in that manipulating the data setincorporates the steps of: a) providing a focus criterion for the imagedata set prior to correction, b) generating a trial artefactcorresponding to a possible artefact cause, c) producing a manipulatedversion of the image data set in which the trial artefact iscounteracted and providing a focus criterion for this version, d)varying the trial artefact and repeating the preceding step iterativelyto produce a manipulated version of the image data set compensated for aversion of the trial artefact associated with optimisation of the focuscriterion to the extent obtainable by iteration, and e) producingcorrected image data from a manipulated version of the image data setassociated with optimisation of the focus criterion.
 2. A methodaccording to claim 1 characterised in that the said trial artefact is afirst trial artefact, and after step (d), steps (b) to (d) are repeatedfor a second trial artefact.
 3. A method according to claim 1characterised in that the image is a physiological image and theunwanted artefact effects are associated with object motion duringoperation of the scanning means (12).
 4. A method according to claim 3characterised in that the scanning means is a magnetic resonance imaging(MRI) scanning means (12).
 5. A method for generating a corrected imageof an object by: (i) employing scanning means (12) to scan an object andderive an image data set from signals received from the object, and (i)manipulating the data set to reduce unwanted artefact effects associatedwith the object and produce corrected image data to yield anartefact-reduced image, characterised in that manipulating the data setincorporates the steps of: a) providing a focus criterion for the imagedata set prior to correction, b) generating a modelled artefactcorresponding to a possible artefact cause, c) producing a manipulatedversion of the image data set in which the modelled artefact iscounteracted and providing a focus criterion for this version, d)varying the modelled artefact and repeating the preceding stepiteratively to produce a manipulated version of the image data setcompensated for a version of the modelled artefact associated withoptimisation of the focus criterion to the extent obtainable byiteration, and e) producing corrected image data from a manipulatedversion of the image data set associated with optimisation of the focuscriterion, wherein the image is a physiological image and the unwantedartefact effects are associated with object motion during operation ofthe scanning means (12), wherein the scanning means is a magneticresonance imaging (MRI) scanning means (12), wherein the focus criterionis an image entropy criterion.
 6. A method according to claim 5characterised in that the artefact effects arise from objectdisplacement and the step of producing a manipulated version of theimage data set incorporates multiplying a plurality of the receivedsignals by respective phase factors.
 7. A method according to claim 6characterised in that the MRI scanning means (12) defines a data readoutdirection, and to correct for a displacement m_(x) ^(j) of an image dataset point in this direction the phase factor is exp[−K_(x);m_(x) ^(j)],where K_(x)=2/FOV_(x)(u−N_(x)/2), FOV_(x) is the field of view in thereadout direction and u is an index number of the said image data setpoint in a line of such points extending in the readout direction andcontaining a number N_(x) of such points.
 8. A method according to claim6 characterised in that the MRI scanning means (12) defines a phaseencode direction, the image data set is a set of lines of data points,and to correct for a displacement m_(y) ^(j) of a line of index number jin the image data set, the phase factor is exp[−K_(y);m_(j) ^(j)], whereK_(y)=2/FOV_(y)(u−N_(y)/2), FOV_(y) is the field of view in the phaseencode direction and N_(y) is the number of lines in the data set.
 9. Amethod for generating a corrected image of an object by: (i) employingscanning means (12) to scan an object and derive an image data set fromsignals received from the object, and (i) manipulating the data set toreduce unwanted artefact effects associated with the object and producecorrected image data to yield an artefact-reduced image, characterisedin that manipulating the data set incorporates the steps of: a)providing a focus criterion for the image data set prior to correction,b) generating a modelled artefact corresponding to a possible artefactcause, c) producing a manipulated version of the image data set in whichthe modelled artefact is counteracted and providing a focus criterionfor this version, d) varying the modelled artefact and repeating thepreceding step iteratively to produce a manipulated version of the imagedata set compensated for a version of the modelled artefact associatedwith optimisation of the focus criterion to the extent obtainable byiteration, and e) producing corrected image data from a manipulatedversion of the image data set associated with optimisation of the focuscriterion, wherein the focus criterion is an image entropy criteriondetermined from sub-divisions of an image which are identified as beingnon-dark.
 10. A method for generating a corrected image of an object by:(i) employing scanning means (12) to scan an object and derive an imagedata set from signals received from the object, and (i) manipulating thedata set to reduce unwanted artefact effects associated with the objectand produce corrected image data to yield an artefact-reduced image,characterised in that manipulating the data set incorporates the stepsof: a) providing a focus criterion for the image data set prior tocorrection, b) generating a modelled artefact corresponding to apossible artefact cause, c) producing a manipulated version of the imagedata set in which the modelled artefact is counteracted and providing afocus criterion for this version, d) varying the modelled artefact andrepeating the preceding step iteratively to produce a manipulatedversion of the image data set compensated for a version of the modelledartefact associated with optimisation of the focus criterion to theextent obtainable by iteration, and e) producing corrected image datafrom a manipulated version of the image data set associated withoptimisation of the focus criterion, wherein the scanning means is anMRI scanning means (12), characterised in that the modelled artefact isassociated with object motion, the image data set is a set of lines ofdata points and the step of iteratively varying the modelled artefact isperformed initially using a data segment incorporating a larger numberof such lines relative to those of a data segment used later.
 11. Amethod for generating a corrected image of an object by: (i) employingscanning means (12) to scan an object and derive an image data set fromsignals received from the object, and (i) manipulating the data set toreduce unwanted artefact effects associated with the object and producecorrected image data to yield an artefact-reduced image, characterisedin that manipulating the data set incorporates the steps of: a)providing a focus criterion for the image data set prior to correction,b) generating a modelled artefact corresponding to a possible artefactcause, c) producing a manipulated version of the image data set in whichthe modelled artefact is counteracted and providing a focus criterionfor this version, d) varying the modelled artefact and repeating thepreceding step iteratively to produce a manipulated version of the imagedata set compensated for a version of the modelled artefact associatedwith optimisation of the focus criterion to the extent obtainable byiteration, and e) producing corrected image data from a manipulatedversion of the image data set associated with optimisation of the focuscriterion, wherein the scanning means comprises an MRI scanning means(12), characterised in that the modelled artefact is associated withobject motion, and the step of iteratively varying the modelled artefactis performed initially using a smaller number of nodes relative to thoseused later.
 12. A method for generating a corrected image of an objectby: (i) employing scanning means (12) to scan an object and derive animage data set from signals received from the object, and (i)manipulating the data set to reduce unwanted artefact effects associatedwith the object and produce corrected image data to yield anartefact-reduced image, characterised in that manipulating the data setincorporates the steps of: a) providing a focus criterion for the imagedata set prior to correction, b) generating a modelled artefactcorresponding to a possible artefact cause, c) producing a manipulatedversion of the image data set in which the modelled artefact iscounteracted and providing a focus criterion for this version, d)varying the modelled artefact and repeating the preceding stepiteratively to produce a manipulated version of the image data setcompensated for a version of the modelled artefact associated withoptimisation of the focus criterion to the extent obtainable byiteration, and e) producing corrected image data from a manipulatedversion of the image data set associated with optimisation of the focuscriterion, wherein the scanning means comprises an MRI scanning means(12), characterised in that the modelled artefact is associated withobject motion, the MRI scanning means (12) is arranged to use a sub-setof the data set to obtain a sub-image, the step of iteratively varyingthe modelled artefact is performed over a region corresponding to thesub-set, and the sub-set is extended later in such iteration.
 13. Animaging system incorporating: (i) scanning means (12) for scanning anobject and deriving an image data set from signals received from theobject, and (ii) processing means (38) for manipulating the image dataset to reduce unwanted artefact effects associated with the object andproduce corrected image data to yield an artefact-reduced image,characterised in that the processing means (38) is arranged to executethe following functions: a) provision of a focus criterion for the imagedata set prior to correction, b) generation of a trial artefactcorresponding to a possible artefact cause, c) production of amanipulated version of the image data set in which the trial artefact iscounteracted and provision of a further focus criterion for thisversion, d) variation of the trial artefact and repetition of thepreceding function iteratively to produce a manipulated version of theimage data set compensated for a version of the trial artefactassociated with optimisation of the focus criterion to the extentobtainable by iteration, and e) production of corrected image data fromthe manipulated version of the image data set associated withoptimisation of the focus criterion.
 14. An imaging system according toclaim 13 characterised in that the said trial artefact is a first trialartefact and after function (d) the processing means (38) is arranged torepeat functions (b) to (d) for a second trial artefact.
 15. An imagingsystem according to claim 13 characterised in that the image is aphysiological image and the unwanted artefact effects are associatedwith object motion during operation of the scanning means (12).
 16. Animaging system according to claim 15 characterised in that the scanningmeans is an MRI scanning means (12).
 17. An imaging systemincorporating: (i) scanning means (12) for scanning an object andderiving an image data set from signals received from the object, and(ii) processing means (38) for manipulating the image data set to reduceunwanted artefact effects associated with the object and producecorrected image data to yield an artefact-reduced image, characterisedin that the processing means (38) is arranged to execute the followingfunctions: a) provision of a focus criterion for the image data setprior to correction, b) generation of a modelled artefact correspondingto a possible artefact cause, c) production of a manipulated version ofthe image data set in which the modelled artefact is counteracted andprovision of a further focus criterion for this version, d) variation ofthe modelled artefact and repetition of the preceding functioniteratively to produce a manipulated version of the image data setcompensated for a version of the modelled artefact associated withoptimisation of the focus criterion to the extent obtainable byiteration, and e) production of corrected image data from themanipulated version of the image data set associated with optimisationof the focus criterion, wherein the image is a physiological image andthe unwanted artefact effects are associated with object motion duringoperation of the scanning means (12), wherein the scanning means is anMRI scanning means (12), wherein the focus criterion is an image entropycriterion.
 18. An imaging system according to claim 17 characterised inthat the artefact effects are associated with object displacement andthe processing means (38) is arranged to implement the function ofproducing a manipulated version of the data set by multiplying aplurality of the received signals by respective phase factors.
 19. Animaging system according to claim 18 characterised in that the MRIscanning means (12) defines a data readout direction, and to correct fora displacement m_(x) ^(j) of an image data set point in this directionthe phase factor is exp[−K_(x).m_(x) ^(j)], whereK_(x)=2/FOV_(x)(u−N_(x)/2), FOV_(x) is the field of view in the readoutdirection and u is an index number of the said image data set point in aline of such points extending in the readout direction and containing anumber N_(x) of such points.
 20. An imaging system according to claim 17characterised in that the MRI scanning means (12) defines a phase encodedirection, the image data set is a set of lines of data points, and tocorrect for a displacement m_(y) ^(j) of a line of index number j in theimage data set the phase factor is exp[−K_(y).m_(y) ^(j)], whereK_(y)=2/FOV_(y)(u−N_(y)/2), FOV_(y) is the field of view in the phaseencode direction and N_(y) is the number of lines in the data set. 21.An imaging system incorporating: (i) scanning means (12) for scanning anobject and deriving an image data set from signals received from theobject, and (ii) processing means (38) for manipulating the image dataset to reduce unwanted artefact effects associated with the object andproduce corrected image data to yield an artefact-reduced image,characterised in that the processing means (38) is arranged to executethe following functions: a) provision of a focus criterion for the imagedata set prior to correction, b) generation of a modelled artefactcorresponding to a possible artefact cause, c) production of amanipulated version of the image data set in which the modelled artefactis counteracted and provision of a further focus criterion for thisversion, d) variation of the modelled artefact and repetition of thepreceding function iteratively to produce a manipulated version of theimage data set compensated for a version of the modelled artefactassociated with optimisation of the focus criterion to the extentobtainable by iteration, and e) production of corrected image data fromthe manipulated version of the image data set associated withoptimisation of the focus criterion, wherein the focus criterion is animage entropy criterion determined from sub-divisions of an image whichare identified as being non-dark.
 22. An imaging system incorporating:(i) scanning means (12) for scanning an object and deriving an imagedata set from signals received from the object, and (ii) processingmeans (38) for manipulating the image data set to reduce unwantedartefact effects associated with the object and produce corrected imagedata to yield an artefact-reduced image, characterised in that theprocessing means (38) is arranged to execute the following functions: a)provision of a focus criterion for the image data set prior tocorrection, b) generation of a modelled artefact corresponding to apossible artefact cause, c) production of a manipulated version of theimage data set in which the modelled artefact is counteracted andprovision of a further focus criterion for this version, d) variation ofthe modelled artefact and repetition of the preceding functioniteratively to produce a manipulated version of the image data setcompensated for a version of the modelled artefact associated withoptimisation of the focus criterion to the extent obtainable byiteration, and e) production of corrected image data from themanipulated version of the image data set associated with optimisationof the focus criterion, wherein the scanning means is an MRI scanningmeans (12), characterised in that the modelled artefact is associatedwith object motion, the image data set is a set of lines of data pointsand the processing means (38) is arranged to perform iterative variationof the modelled artefact using initially a data segment incorporating alarger number of such lines relative to those of a data segment usedlater in the iteration function.
 23. An imaging system incorporating:(i) scanning means (12) for scanning an object and deriving an imagedata set from signals received from the object, and (ii) processingmeans (38) for manipulating the image data set to reduce unwantedartefact effects associated with the object and produce corrected imagedata to yield an artefact-reduced image, characterised in that theprocessing means (38) is arranged to execute the following functions: a)provision of a focus criterion for the image data set prior tocorrection, b) generation of a modelled artefact corresponding to apossible artefact cause, c) production of a manipulated version of theimage data set in which the modelled artefact is counteracted andprovision of a further focus criterion for this version, d) variation ofthe modelled artefact and repetition of the preceding functioniteratively to produce a manipulated version of the image data setcompensated for a version of the modelled artefact associated withoptimisation of the focus criterion to the extent obtainable byiteration, and e) production of corrected image data from themanipulated version of the image data set associated with optimisationof the focus criterion, wherein the scanning means is an MRI scanningmeans (12), characterised in that the modelled artefact is associatedwith object motion and the processing means (38) is arranged to performiterative variation of the modelled artefact using initially a smallernumber of nodes relative to those used later in the iteration function.24. An imaging system incorporating: (i) scanning means (12) forscanning an object and deriving an image data set from signals receivedfrom the object, and (ii) processing means (38) for manipulating theimage data set to reduce unwanted artefact effects associated with theobject and produce corrected image data to yield an artefact-reducedimage, characterised in that the processing means (38) is arranged toexecute the following functions: a) provision of a focus criterion forthe image data set prior to correction, b) generation of a modelledartefact corresponding to a possible artefact cause, c) production of amanipulated version of the image data set in which the modelled artefactis counteracted and provision of a further focus criterion for thisversion, d) variation of the modelled artefact and repetition of thepreceding function iteratively to produce a manipulated version of theimage data set compensated for a version of the modelled artefactassociated with optimisation of the focus criterion to the extentobtainable by iteration, and e) production of corrected image data fromthe manipulated version of the image data set associated withoptimisation of the focus criterion, wherein the scanning means is anMRI scanning means (12), characterised in that the modelled artefact isassociated with object motion, the MRI scanning means (12) is arrangedto use a sub-set of the data set to obtain a sub-image, the processingmeans (38) is arranged to perform iterative variation of the modelledartefact over a region corresponding to the sub-set and to extend thesub-set later in the iteration function.